Impacts of climate change on the potential distribution of Pulex simulans and Polygenis gwyni

Abstract Pulex simulans and Polygenis gwyni are vectors of many flea‐borne diseases. They were widely recorded in the United States and Mexico between 1970 and 2000. Maximum entropy models were used to explore the habitats of both fleas under different climate scenarios to provide the scientific basis for the surveillance and control of flea‐borne diseases. We screened climate variables by principal component analysis and Pearson's correlation test and evaluated model performance by ROC curve. ArcMap was used to visualize expressions. Under current climatic conditions, the medium and highly suitable areas for P. simulans are estimated to be 9.16 × 106 km2 and 4.97 × 106 km2, respectively. These regions are predominantly located in South America, along the Mediterranean coast of Europe, the southern part of the African continent, the Middle East, North China, and Australia. For P. gwyni, the medium and highly suitable areas under current climatic conditions are approximately 4.01 × 106 and 2.04 × 106 km2, respectively, with the primary distribution in North China extending to the Himalayas, near the Equator in Africa, and in a few areas of Europe. Under future climate scenarios, in the SSP3‐7.0 scenario for the years 2081–2100, the area of high suitability for P. simulans is projected to reach its maximum. Similarly, in the SSP2‐4.5 scenario for 2061–2080, the area of high suitability for P. gwyni is expected to reach its maximum. Under global climate change, there is a large range in the potential distribution for both fleas, with an overall upward trend in the area of habitat under future climate scenarios. Governments should develop scientific prevention and control measures to prevent the invasive alien species flea.


| INTRODUC TI ON
The advancement of human economic activities has been a significant driver of climate change (Trenberth, 2018), which in turn has heightened the risk of infectious disease transmission.The escalation in global temperatures has broadened the geographical distribution of arthropod vectors, enabling their migration to higher latitudes and altitudes, thereby intensifying the dissemination of vectorborne diseases (Bale et al., 2002;Githeko et al., 2000).Among these vectors, P. simulans and P. gwyni are noteworthy for their role in the transmission of plague.P. simulans is thought to be able to carry the Yersinia pestis, and P. gwyni, a flea species associated with cotton rats, has facilitated the spread of plague to Sigmodon hispidus and Rattus norveicus in the southern and southwestern United States (Holdenried, 1952;Poje et al., 2020).
In terms of distribution, research by AG Mutlow has highlighted that P. simulans exhibits a broad geographic presence across the United States, Central, and South America, with a diverse host range (Mutlow et al., 2006).Similarly, Durden, LA et al.'s study identified the primary habitat of P. gwyni in the Carolinas along the southeastern Atlantic coast of the United States (Durden et al., 1999), underscoring the ecological and epidemiological significance of these flea species in the context of disease transmission.Climate factors are closely related to the occurrence of vector-borne infectious diseases.Climate change greatly affects the behavior of flea species and their ability to transmit diseases (Gage et al., 2008).For example, factors such as temperature, precipitation, and relative humidity can affect the growth and survival of fleas (Guernier et al., 2014), and research by Xu and his team on the plague in China during the 19th and 20th centuries found a nonlinear relationship between the intensity of plague and precipitation (Xu et al., 2011).Climate factors also affect the population numbers and distribution of vectors (Liu et al., 2018), creating new epidemic risks.
Therefore, using ecological niche models, which are particularly suitable for ecology and biogeography, to assess the environmental conditions within a species' distribution range helps predict the possibility of biological invasions and plays an important role in preventing outbreaks of flea-borne diseases.The Maximum Entropy (MaxEnt) model (Feng et al., 2019;Peterson, 2003) stands out in this domain.MaxEnt, known for its continuous prediction capabilities (Phillips et al., 2006), is adept at generating various predictive metrics such as omission rates, receiver operating characteristic (ROC) curves, and response curves.This model is acclaimed for its exceptional predictive accuracy, as highlighted by Elith et al. (2006), making it a preferred choice in geographic modeling.Further emphasizing MaxEnt's capabilities, one research (Zhang & Wang, 2022) demonstrates its superiority over other models like logistic regression, artificial neural networks, and support vector machines.The MaxEnt model's strengths lie in its ability to account for interactions and the significance of spatial variables, coupled with generating comprehensible response curves, thereby offering a more coherent explanation of ecological phenomena.This robust framework renders MaxEnt an invaluable tool for understanding and predicting the dynamics of species distribution, particularly in the context of changing climatic conditions.
As economic globalization, urbanization, and increases in personal income fuel trade between nations and travel among individuals (Wu et al., 2017), along with climate change, there has been a facilitated widespread dissemination of pathogens across the globe.Consequently, there is an elevated risk of diseases transmitted by fleas.The innovation of this study lies in its exploration of the changes in suitable habitats for the plague bacterium-carrying vector species P. simulans and P. gwyni, driven by current and future climate scenarios, through the use of the Maxent model.Although these two fleas have not been found in other areas for the time being, they are at risk of invasion under certain conditions.So, this investigation is crucial for preventing the invasion of these vectors and, by extension, is vital for the prevention of flea-borne diseases, holding significant public health implications.

| Distribution point data acquisition
The data of fleas were downloaded from Global Biodiversity Information Facility (https:// www.gbif.org; accessed on 23 March 2023).We checked the accuracy of the data by selecting only distribution point data from 1970 to 2000 and deleting blank and duplicate data.In ENMTools software (https:// www.activ estate.com/ produ cts/ perl/ downl oads; accessed on 15 April 2023), click Trim duplicate occurrences, import the csv file with distribution points, and select Grid cell, then we can get the csv file with deleted redundant data in the same grid.Finally, 50 points for P. simulans and 17 points for P. gwyni were obtained and then stored as a CSV file containing species names, longitude, and latitude.The distribution can be seen in Figure 1.  1) and selected the environmental variables for P. gwyni through a correlation significance test (Table 2).Variables with an absolute correlation coefficient >.8 were excluded if they contributed less.Furthermore, to ascertain the presence of multicollinearity, the variance inflation factor (VIF) was calculated for all variables, with all VIF values being below 10.

| Parameter optimization for maximum entropy models
MaxEnt uses a maximum entropy model to predict the likely range of species from known species distribution information and environmental data.However, the predictions of unoptimized models can be subject to significant fit bias, leading to a mis-assessment of species' ecological niches (Kong et al., 2019).Therefore, the parameters need to be changed to optimize the model.The regularization multiplier has eight levels: 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4. Feature classes have five parameters including Linear (L), Quadratic (Q), Product (P), Threshold (T), Hinge (H)."Import to Gridfile" with the type set to BIL, the BIL format files were converted into grd and gil formats.The best model was determined using the Akaike Information Criterion (AIC) in the ENMeval package in R software, where an optimal model is indicated by a deltaAICc of 0. In this study, the optimal model for P. simulans was identified with a regularization multiplier (rm) of 1 and feature class (fc) as QHP; for P. gwyni, the optimal model had a rm of 2 and fc as LQH.The training and testing datasets are randomly generated, with 75% of the data points used for training and 25% used for testing.Cross-validation with the testing and validation sets is utilized to avoid overfitting.

| Distribution and change of suitable area
The suitable area, which refers to the range most conducive to the growth and development of an organism, for P. simulans and P. gwyni  2) was used to test the accuracy of results, and the AUC value >0.9 indicates high accuracy.

| Distribution concerning environmental variables
The ROC curve (Figure 2) is used to test the accuracy of results, and the AUC value greater than 0.9 indicates high accuracy and the Jackknife test (Figure 3) is used to examine the contribution of each variable to influencing distribution (Çoban et al., 2020).We predicted the distribution of potential suitable area for P. simulans and P. gwyni under current and future climatic conditions by parameter-optimized maxent.In our study, we identified 13 environmental variables contributing more than 2% to the distribution of P. simulans and performed a PCA to synthesize three new variables.Meanwhile, we selected five environmental variables that contributed more than 2% to the distribution of P. gwyni, each with a correlation coefficient lower than 0.8.
Their VIF values were below 10.These environmental variables were incorporated into a parameter-optimized MaxEnt model to predict their potential suitable habitat distributions under current and future climate scenarios.In the modeling, it was observed that the key determinants of P. simulans habitat are primarily PCA1 (represented mainly by precipitation in August and September) and PCA3 (represented mainly by the highest temperatures in February and July), whereas the crucial factors for P. gwyni habitat were identified as September precipitation and the minimum temperature in January.

| Distribution of potentially suitable areas under current climate scenarios
According to the results derived from MaxEnt, as illustrated in | 5 of 10 WANG et al.

| Distribution of potentially suitable areas under future climate scenarios
The findings indicate that among 16 future climate scenarios, the total suitable habitat area for P. simulans increased under six scenarios, whereas for P. gwyni, an increase was observed under eight scenarios.See Table 3 for all areas.
As shown in Figure 6, the maximum area of high suitability for

| DISCUSS ION
Fleas are vectors of some animal-borne diseases.As the global climate changes, which significantly affects the growth and habitat range of fleas and is more likely to transmit diseases, fleas can transmit pruritic papules of the ankle, rabbit fever, cat-scratch disease, bartonellosis, rickettsiosis, and plague (Bitam et al., 2010) and have become a problem of public health importance.We find through the model that temperature and precipitation play important roles in the distribution of these two flea species.As seen by Kreppel et al. (2016), temperature notably affects the development time of flea larvae and pupae.Fleas are highly sensitive to temperature changes, particularly during their developmental stages (van der Mescht et al., 2016), where lower temperatures in winter may reduce flea survival rates (Blagburn & Dryden, 2009;Silverman et al., 1981).
Environmental stability can enhance the stability of flea populations, thereby increasing their numbers and activity.This offers more opportunities for contact with hosts, elevating the risk of biological invasions.Precipitation's impact on flea populations is equally critical; anomalous abundance of rainfall can lead to increased humidity, fostering the proliferation of fungi and mites within rodent burrows, ultimately leading to a decrease in flea populations (Krasnov, 2008;Wimsatt & Biggins, 2009).and Asia respond most significantly to temperature differences (Krasnov et al., 2020).
ENMeval is the first R package to optimize the quantitative complexity assessment for MaxEnt (Kass et al., 2021).What stands out about our research lies in the use of ENMeval in R to optimize the  et al. (2014) find that among land use types, the flea index is higher in fallow and natural forests, while the lowest flea index is found in plantation monocultures and annual mixed crops, and does not incorporate the effect of flea hosts on their abundance (Russell et al., 2018), BR Krasnov finds that the similarity between habitats of fleas composition increases with increasing similarity in hosts composition (Krasnov et al., 2006), all of which can potentially influence flea distribution and hence predictions.Despite the absence of host-related data, we believe that the predicted distributions of P. simulans and P. gwyni are relatively consistent with their host preferences.Given the broad range of hosts for P. simulans (Mutlow et al., 2006), there is a potential risk of invasion; the Rattus norvegicus, one of the hosts for P. gwyni, is widely distributed globally (Worth, 1950), thus posing a potential invasion risk as well, As it is difficult to distinguish morphologically between female P. simulans and P. irritans, it may be possible to misidentify flea species to some extent, which may affect the accuracy of our predictions, we only download data in WorldClim from 1970 to 2000 under current climate conditions, which may have affected predictions due to data limitations.I believe that future research should not only consider climatic factors but also other environmental variables, and must take into account the parasitic relationships of organisms to predict the distribution of their suitable habitats.

| CON CLUS IONS
We used the MaxEnt to predict the potential habitats of P. simulans and P. gwyni and found that under current climatic conditions, the Additionally, stringent border health checks should be enforced in other suitable regions to prevent the introduction of vector species.This approach is vital for mitigating the health risks associated with climate-induced changes in the distribution of flea vectors.
Environmental data from 1970 to 2000 were downloaded from WorldClim (version 2.1, https:// world clim.org/ ; accessed on: 24 March 2023), including 19 Bioclimatic variables, 36 historical monthly climate data and elevation data, with a spatial resolution of 5 mm.and download environmental data for 2021-2040, 2041-2060, 2061-2080, 2081-2100 for the four Shared Socio-economic Pathways (SSPs): 126, 245, 370 and 585 scenarios.All environmental data were downloaded in TIFF format and then converted to ASCII format using ArcGIS (version 10.5), purchased by the Chinese Center for Disease Control and Prevention.The filtered distributions of P. simulans and P. gwyni and the 56 environmental variables were imported in MaxEnt software (version 3.4.1,https:// biodi versi tyinf ormat ics.amnh.org; accessed on 13 April 2023) and ran once to obtain the variable contributions.These environmental variables with a contribution rate greater than or equal to | 3 of 10 WANG et al. 2% were screened out.Statistical analyses revealed high correlations among certain environmental variables at distribution points.To clarify the similarities among bioclimatic variables, we initially conducted a Kaiser-Meyer-Olkin (KMO) test.A Kaiser-Meyer-Olkin Measure of Sampling Adequacy greater than 0.6 indicates the suitability of employing the principal component analysis (PCA).The environmental variables of the distribution points for P. simulans yielded a Kaiser-Meyer-Olkin Measure of Sampling Adequacy of 0.669, whereas those for P. gwyni recorded a value of 0.458.And Bartlett's test of sphericity yields a significance level of <.001.Consequently, we utilized PCA to process the environmental variables at the P. simulans distribution points (Table

R
software and DIVA-GIS correlated the environmental data and adjusted the parameters to optimize the MaxEnt model.Environmental data were correlated using ArcGIS, DIVA-GIS, and R software.Initially, raster files were exported from ArcGIS in the BIL format.In DIVA-GIS, by selecting "Data" from the menu and choosing F I G U R E 1 Occurrence points of Pulex simulans and Polygenis gwyni.TA B L E 1 Results of principal component analysis.

Figure 4 ,
Figure4, under current climatic conditions, the potentially suitable habitat area for species P. simualns spans 36.32 × 10 6 km 2 .The areas of medium and high suitability are respectively 9.16 × 10 6 and 4.97 × 10 6 km 2 .These regions are predominantly located in South America, including Peru, Argentina, Chile, and Brazil; in Europe along the Mediterranean coast; in Africa, particularly in the southern parts of the continent such as South Africa and Madagascar; and in other areas of medium to high suitability encompass the Middle East, North China, and Australia.As depicted in Figure5, under the current climate conditions, the potential suitable habitat area for species P. gwyni is 15.67 × 10 6 km 2 , with the medium and high suitability areas measuring 4.01 × 10 6 and 2.04 × 10 6 km 2 , respectively.These regions are mainly distributed from the North China region to the Himalayan range in Asia, near the equator in Africa, and in a few parts of Europe.

P
. simulans across different time periods is observed in the SSP3-7.0scenario for the years 2081-2100, where the high suitability area reached its peak at 5.64 × 10 6 km 2 , compared to the currently suitable area of 4.97 × 10 6 km 2 .In contrast, under the SSP5-8.5 scenario for the years 2081-2100, the high suitability area is at its minimum, measuring 4.71 × 10 6 km 2 .F I G U R E 3 The results of the Jackknife test of variable importance for Pulex simulans (a) and Polygenis gwyni (b): The green bands indicate the degree of specificity of the variable, with shorter bars indicating the variable contains unique information and is more likely to influence the distribution of the species.The blue bands indicate the variable's effectiveness on the species' distribution, with longer bars indicating that the variable contains more effective information.All bioclimatic variables utilized in the model are deemed essential.F I G U R E 4 Suitable areas for Pulex simulans under current climate scenarios.For P. gwyni, as illustrated in Figure7, the largest area of high suitability across different time periods is noted in the SSP2-4.5 scenario for the years 2061-2080, where the high suitability area is the highest at 2.73 × 10 6 km 2 , in comparison to the current suitability area of 2.04 × 10 6 km 2 .The smallest area of high suitability is observed in the year 2061 under the SSP5-8.5 scenario, measuring 1.04 × 10 6 km 2 .
Results obtained through Maxent clearly indicate significant changes in the distribution regions of P. simulans and P. gwyni, with such changes also dependent on various climate scenarios.Compared to P. gwyni, P. simulans has a larger area of current and future habitat; the potential habitat of P. simulans is South America; Europe along the Mediterranean coast; Africa, particularly in the southern parts of the continent; other areas of medium to high suitability encompass the Middle East, North China, and Australia.Differences in climate change intensity significantly influence the global distribution of their habitats, showcasing the most notable expansion during 2061-2080 under the SSP1-2.6 scenario.The potential suitable areas for P. gwyni mainly encompass the North China region to the Himalayan range in Asia, regions near the equator in Africa, and several parts of Europe.Previous studies (Li et al., 2023) have also explored the impact of climate change on the potential suitable habitats of fleas, with similarities to our study in that we both employed the Maxent model.Their research predicted the distribution of two flea species found in North America, Peromyscopsylla hesperomys and Orchopeas sexdentatus, which have distribution areas broadly similar to ours, including South America, the Mediterranean coast, and the North China region.The key difference lies in the primary environmental factors affecting the distribution of fleas; precipitation had a relatively minor impact on the flea species studied by them compared to those in our research.Earlier research by Boris R. Krasnov found that fleas in North America show the strongest response to variations in precipitation, whereas fleas in South America, Europe,

F
I G U R E 5 Suitable areas for Polygenis gwyni under current climate scenarios.parameters and avoid over-fitting of the ecological niche model (ENM), resulting in more accurate results, and in the fact that the study predicts the distribution of the two fleas from a global perspective, which has broader public health implications.The weaknesses are that the ecological data do not include, e.g., normalized vegetation index, land use, relationship with hosts, etc. Young et al. (2015) find that vegetation cover negatively affects the intensity of flea infestation on hosts, which affects the intensity of flea infestation of hosts; Hieronimo

F
I G U R E 6 Suitable areas for Pulex simulans under future climate scenarios (a) SSP245:2021-2040; (b) SSP585:2041-2060; (c) SSP126:2061-2080; (d) SSP370:2081-2100.F I G U R E 7 Suitable areas for Polygenis gwyni under future climate scenarios (a) SSP126:2021-2040; (b) SSP585:2041-2060; (c) SSP245:2061-2080; (d) SSP370:2081-2100.medium and high suitability zones for P. simualns are mainly located in South America, the Mediterranean coast of Europe, the southern part of the African continent, the Middle East, North China and Australia, while the medium and high suitability zones for P. gwyni are mainly located in Asia, from North China to the Himalayan mountain ranges, in Africa near the equator, and in a few parts of Europe.Climate change is expanding the distribution ranges of P. simulans and P. gwyni, making it imperative to proactively address the increased risk of flea-borne diseases driven by these changes.Enhanced surveillance in areas highly suitable for the survival of these vectors is essential to prevent disease outbreaks.
Global habitat distribution of Pulex simulans and Polygenis gwyni under current and future climate conditions (×10 6 km 2 ).